#############################Appendix Analyses########################

#Last run by authors on macOS Monterey v 12.1
# Rstudio verison 2022.02.3 
# Running R version 4.1.2 Bird Hippie

#Define functions
#Produce figures from Appendix

#Files In: data_panel.csv, engagement_panel.csv

###########################################################################

###Required packages
#install.packages('pacman')
pacman::p_load(plm, lmtest, sandwich, tidyverse, tidyr, broom, stargazer,
               fixest, this.path, bacondecomp, did, ggplot2, forcats,
               dplyr, plyr, lubridate, knitr, kableExtra)

#sets working directory to current folder
setwd(dirname(this.path()))

#pulls required data and functions (both for main analysis and for appendix)
source('VZ_source_code.R')


################################
###Figure A1: Expert Validation

prop_table<-read.csv('ReplicationData/prop_table.csv')

ggplot(prop_table, aes(fill=Rating, y=Proportion, x=Topic))+ 
  geom_bar(position="dodge", stat="identity")+scale_fill_grey()+
  theme_minimal(base_size=18)
ggsave("plots/tweet_topic_validation.pdf", width = 11, height = 8)

######################################
#Figures A2: Weekly volume of tweets

weektweets<-read.csv('ReplicationData/weektweets.csv')
weektweets$week<-as.Date(weektweets$week)

ggplot(subset(weektweets, outcome=='foreign policy'), aes(x=week, y=tpd)) + 
  geom_line() + 
  labs(y = "Weekly Volume of Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/foreign_policy_tweets_weekly_combined.pdf", width = 11, height = 7)

ggplot(subset(weektweets, outcome=='service'), aes(x=week, y=tpd)) + 
  geom_line() + 
  labs(y = "Weekly Volume of Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/service_tweets_weekly_combined.pdf", width = 11, height = 7)


ggplot(subset(weektweets, outcome=='protest'), aes(x=week, y=tpd)) + 
  geom_line() + 
  labs(y = "Weekly Volume of Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/protest_tweets_weekly_combined.pdf", width = 11, height = 7)

ggplot(subset(weektweets, outcome=='criticism'), aes(x=week, y=tpd)) + 
  geom_line() + 
  labs(y = "Weekly Volume of Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/harsh_criticism_tweets_weekly_combined.pdf", width = 11, height = 7)


######################################
#Figures A3: Mean tweets


subset(weektweets, outcome=='foreign policy') %>% 
  ggplot() + 
  aes(x = week, y = mean_tpd) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/foreign_policy_mean_tweets_weekly_combined.pdf", width = 11, height = 7)


subset(weektweets, outcome=='service') %>% 
  ggplot() + 
  aes(x = week, y = mean_tpd) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/service_mean_tweets_weekly_combined.pdf", width = 11, height = 7)


subset(weektweets, outcome=='protest') %>% 
  ggplot() + 
  aes(x = week, y = mean_tpd) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/protest_mean_tweets_weekly_combined.pdf", width = 11, height = 7)


subset(weektweets, outcome=='criticism') %>% 
  ggplot() + 
  aes(x = week, y = mean_tpd) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/harsh_criticism_mean_tweets_weekly_combined.pdf", width = 11, height = 7)



######################################
#Figures A4: Mean tweets


subset(weektweets, outcome=='foreign policy') %>% 
  ggplot() + 
  aes(x = week, y = prop) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/foreign_policy_prop_tweets_weekly_combined.pdf", width = 11, height = 7)


subset(weektweets, outcome=='service') %>% 
  ggplot() + 
  aes(x = week, y = prop) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/service_prop_tweets_weekly_combined.pdf", width = 11, height = 7)


subset(weektweets, outcome=='protest') %>% 
  ggplot() + 
  aes(x = week, y = prop) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/protest_prop_tweets_weekly_combined.pdf", width = 11, height = 7)


subset(weektweets, outcome=='criticism') %>% 
  ggplot() + 
  aes(x = week, y = prop) + 
  geom_line() + 
  labs(y = "Mean Weekly Tweets", x = "Date")+
  theme_minimal(base_size=22)
ggsave("plots/harsh_criticism_prop_tweets_weekly_combined.pdf", width = 11, height = 7)



####################################
#####Figures A5 A6: Classified Tweets

relevant_military<-getting_models('perc_relevant_military', pdf)
positive_military<-getting_models('perc_positive_military', pdf)
negative_military<-getting_models('perc_negative_military', pdf)

relevant_sanctions<-getting_models('perc_relevant_sanctions', pdf)
positive_sanctions<-getting_models('perc_positive_sanctions', pdf)
negative_sanctions<-getting_models('perc_negative_sanctions', pdf)

models<-list(relevant_military, positive_military, negative_military)

dv_names<-c('Relevant', 'Positive', 'Negative')

dv_types<-c('bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

military_f8_coefplot<-coef_plot(models, dv_names, dv_types, model_types, 'military_f8', 'Classified Intervention Tweets and Exile')



models<-list(relevant_sanctions, positive_sanctions, negative_sanctions)

dv_names<-c('Relevant', 'Positive', 'Negative')

dv_types<-c('bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'sanctions_f8', 'Classified Sanctions and Exile')


##############################################
###############Figure A7: References to exile


#Event study
es_exile <- plm(perc_exile ~ lead_lags+month+log(num_tweets+1), data = pdf, model = "within")
es_exile_coef <- coeftest(es_exile, vcov = vcovHC(es_exile, type = "HC1", cluster = "group"))
es_exile_plot<-es_plot(es_exile_coef, 'Change in % Exile Tweets')
ggsave("plots/es_exile_plot.pdf", es_exile_plot, width = 8, height = 7)


#Two way fixed effects
exile_models<-getting_models('perc_exile', pdf)

models<-exile_models

dv_names<-c('Exile')

dv_types<-c('bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'exile', 'Exile References and Exile')


#####################################
##########Figure A8: results by location of exile

md_cont_fp <- plm(perc_foreign_policy~ tweeted_exile*US_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_cont_fp <- coeftest(md_cont_fp, vcov = vcovHC(md_cont_fp, type = "HC1", cluster = "group"))

md_cont_sv <- plm(perc_service~ tweeted_exile*US_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_cont_sv <- coeftest(md_cont_sv, vcov = vcovHC(md_cont_sv, type = "HC1", cluster = "group"))

md_cont_pr <- plm(perc_protest~ tweeted_exile*US_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_cont_pr <- coeftest(md_cont_pr, vcov = vcovHC(md_cont_pr, type = "HC1", cluster = "group"))

md_cont_cr <- plm(perc_harsh_criticism~ tweeted_exile*US_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_cont_cr <- coeftest(md_cont_cr, vcov = vcovHC(md_cont_cr, type = "HC1", cluster = "group"))

###colombia

md_co_cont_fp <- plm(perc_foreign_policy~ tweeted_exile*CO_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_co_cont_fp <- coeftest(md_co_cont_fp, vcov = vcovHC(md_co_cont_fp, type = "HC1", cluster = "group"))

md_co_cont_sv <- plm(perc_service~ tweeted_exile*CO_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_co_cont_sv <- coeftest(md_co_cont_sv, vcov = vcovHC(md_co_cont_sv, type = "HC1", cluster = "group"))

md_co_cont_pr <- plm(perc_protest~ tweeted_exile*CO_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_co_cont_pr <- coeftest(md_co_cont_pr, vcov = vcovHC(md_co_cont_pr, type = "HC1", cluster = "group"))

md_co_cont_cr <- plm(perc_harsh_criticism~ tweeted_exile*CO_destination+month+log(num_tweets+1), data = pdf, model = "within")
coef_md_co_cont_cr <- coeftest(md_co_cont_cr, vcov = vcovHC(md_co_cont_cr, type = "HC1", cluster = "group"))

coef_df_fp<-subset(data.frame(tidy(coef_md_cont_fp), coef.int=T), term=='tweeted_exile:US_destinationus')
coef_df_sv<-subset(data.frame(tidy(coef_md_cont_sv), coef.int=T), term=='tweeted_exile:US_destinationus')
coef_df_pr<-subset(data.frame(tidy(coef_md_cont_pr), coef.int=T), term=='tweeted_exile:US_destinationus')
coef_df_cr<-subset(data.frame(tidy(coef_md_cont_cr), coef.int=T), term=='tweeted_exile:US_destinationus')

coef_df_cofp<-subset(data.frame(tidy(coef_md_co_cont_fp), coef.int=T), term=='tweeted_exile:CO_destinationcolombia')
coef_df_cosv<-subset(data.frame(tidy(coef_md_co_cont_sv), coef.int=T), term=='tweeted_exile:CO_destinationcolombia')
coef_df_copr<-subset(data.frame(tidy(coef_md_co_cont_pr), coef.int=T), term=='tweeted_exile:CO_destinationcolombia')
coef_df_cocr<-subset(data.frame(tidy(coef_md_co_cont_cr), coef.int=T), term=='tweeted_exile:CO_destinationcolombia')

coef_df<-rbind(coef_df_fp, coef_df_sv, coef_df_pr, coef_df_cr,
               coef_df_cofp, coef_df_cosv, coef_df_copr, coef_df_cocr)

coef_df$conf.low<-coef_df$estimate-coef_df$std.error*1.96
coef_df$conf.high<-coef_df$estimate+coef_df$std.error*1.96
coef_df$model<-c(rep('Exiled in US', 4), rep('Exiled in CO', 4))
coef_df$name<-rep(c('Foreign Policy', 'Service Provision',
                    'Protest', 'Criticism'), 2)
coef_df$name<-rev(factor(coef_df$name, levels=as.character(unique(coef_df$name))))


bydestination_effects_plot<-ggplot(coef_df, aes(shape=model,linetype=model, x=name, y=estimate))+
  geom_pointrange(aes(ymin=conf.low, ymax=conf.high), position=position_dodge(width=.3))+
  theme_bw()+
  theme(legend.position='bottom', axis.title=element_text(size=16),
        text=element_text(size=18))+
  coord_flip()+
  geom_hline(aes(yintercept=0), lty=4)+
  scale_y_continuous(name='Change in % tweets, relative to other exiles')+
  scale_shape_discrete(name=' ', guide=guide_legend(reverse=T))+
  scale_linetype_discrete(name=' ', guide=guide_legend(reverse=T))+
  xlab("")+scale_x_discrete(labels=unique(coef_df$name))+
  theme(plot.title = element_text(hjust = 0.5, size=18))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
ggsave("plots/bydestination_effects_plot.pdf", bydestination_effects_plot, width = 8.5, height = 6.5)



sgd<-stargazer(list(coef_md_cont_fp, coef_md_cont_sv, coef_md_cont_pr, coef_md_cont_cr),
               omit=c('month', 'trend', 'actor.id', 'Constant'),
               title='Effects by Destination',
               align=TRUE, header=FALSE, label='table:bydestination',
               column.labels = dv_names,
               covariate.labels=c('Tweeted from Exile', 'Number of Tweets', 
                                  'Tweeted from Exile*US Destination'),
               omit.stat=c('f', 'adj.rsq'),
               model.numbers=F,
               add.lines=list(c("Month FEs", 'Y', 'Y', 'Y', "Y"),
                              c('User FEs', 'Y', 'Y', 'Y', 'Y')),
               font.size = 'footnotesize',
               column.sep.width = "-10pt",
               dep.var.labels.include=FALSE, dep.var.caption = "", no.space=TRUE,
               notes=c("Robust standard errors clustered at the user level."))

sgd2<-stargazer(list(coef_md_co_cont_fp, coef_md_co_cont_sv, coef_md_co_cont_pr, coef_md_co_cont_cr),
                omit=c('month', 'trend', 'actor.id', 'Constant'),
                title='Effects by Destination',
                align=TRUE, header=FALSE, label='table:bydestination',
                column.labels = dv_names,
                covariate.labels=c('Tweeted from Exile', 'Number of Tweets', 
                                   'Tweeted from Exile*CO Destination'),
                omit.stat=c('f', 'adj.rsq'),
                model.numbers=F,
                add.lines=list(c("Month FEs", 'Y', 'Y', 'Y', "Y"),
                               c('User FEs', 'Y', 'Y', 'Y', 'Y')),
                
                font.size = 'footnotesize',
                column.sep.width = "-10pt",
                dep.var.labels.include=FALSE, dep.var.caption = "", no.space=TRUE,
                notes=c("Robust standard errors clustered at the user level."))


write(c(sgd[1:(length(sgd)-5)], '\\end{tabular} ', sgd2[5:length(sgd2)]), paste('tables/', 'bydestination_effects', '_table.tex', sep=''))





#########################################
###########Figure A9: T-test results

ttest_results<-read.csv('ReplicationData/ttest_results.csv')
ttest_results$period<-factor(ttest_results$period)

ttest_result_plot<-ggplot(ttest_results, aes(shape=period,linetype=period, x=label, y=est))+
  geom_pointrange(aes(ymin=conf.low, ymax=conf.high), position=position_dodge(width=.6), size=.8)+
  theme_bw()+theme(legend.position='bottom', axis.title=element_text(size=22),
                   text=element_text(size=22),
                   axis.text=element_text(size=22))+
  coord_flip()+ggtitle('')+geom_hline(aes(yintercept=0), lty=4)+
  scale_y_continuous(name='Difference in % Tweets')+
  scale_shape_discrete(name='Window (in Days)', guide=guide_legend(reverse=T))+
  scale_linetype_discrete(name='Window (in Days)', guide=guide_legend(reverse=T))+
  xlab("")+
  theme(plot.title = element_text(hjust = 0.5, size=18))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))

ggsave("plots/ttest_result_plot.pdf", ttest_result_plot, width = 10, height = 10)



#############################################
########Figure B1: Including months with no tweets

foreign_policy_wzeros<-getting_models('perc_foreign_policy', pdf_wzeros)
protest_wzeros<-getting_models('perc_protest', pdf_wzeros)
criticism_wzeros<-getting_models('perc_harsh_criticism', pdf_wzeros)
services_wzeros<-getting_models('perc_service', pdf_wzeros)


models<-list(foreign_policy_wzeros, services_wzeros, protest_wzeros, criticism_wzeros)

dv_names<-c('Foreign Policy','Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'wzeros', 'Including Months With No Tweets')




############################################
######Figures B2-B5: Alternative Event Study
pdf$lead_lags_b<-factor(pdf$lead_lags_b, levels=c('month_before', setdiff(unique(pdf$lead_lags_b), 'month_before')))

#foreign policy
es_foreign_policy_b <- plm(perc_foreign_policy ~ lead_lags_b+month+log(num_tweets+1), data = pdf, model = "within")
es_foreign_policy_b_coef <- coeftest(es_foreign_policy_b, vcov = vcovHC(es_foreign_policy_b, type = "HC1", cluster = "group"))
es_foreign_policy_b_plot<-es_plot(es_foreign_policy_b_coef, 'Change in % Foreign Action Tweets')
ggsave("plots/es_foreign_policy_b_plot.pdf", es_foreign_policy_b_plot, width = 8.5, height = 6.5)


#protest
es_protest_b <- plm(perc_protest ~ lead_lags_b+month+log(num_tweets+1), data = pdf, model = "within")
es_protest_b_coef <- coeftest(es_protest_b, vcov = vcovHC(es_protest_b, type = "HC1", cluster = "group"))
es_protest_b_plot<-es_plot(es_protest_b_coef, 'Change in % Protest Tweets')
ggsave("plots/es_protest_b_plot.pdf", es_protest_b_plot, width = 8.5, height = 6.5)


#foreign policy
es_criticism_b <- plm(perc_harsh_criticism ~ lead_lags_b+month+log(num_tweets+1), data = pdf, model = "within")
es_criticism_b_coef <- coeftest(es_criticism_b, vcov = vcovHC(es_criticism_b, type = "HC1", cluster = "group"))
es_criticism_b_plot<-es_plot(es_criticism_b_coef, 'Change in % Tweets on Harsh Criticism')
ggsave("plots/es_criticism_b_plot.pdf", es_criticism_b_plot, width = 8, height = 7)

#services
es_services_b <- plm(perc_service ~ lead_lags_b+month+log(num_tweets+1), data = pdf, model = "within")
es_services_b_coef <- coeftest(es_services_b, vcov = vcovHC(es_services_b, type = "HC1", cluster = "group"))
es_services_b_plot<-es_plot(es_services_b_coef, 'Change in % Tweets on Service Provision')
ggsave("plots/es_services_b_coef.pdf", es_services_b_plot, width = 8, height = 7)


#############################################
############Figure B6: Including only politicians

foreign_policy_pols<-getting_models('perc_foreign_policy', subset(pdf, type=='politician'))
protest_pols<-getting_models('perc_protest', subset(pdf, type=='politician'))
criticism_pols<-getting_models('perc_harsh_criticism', subset(pdf, type=='politician'))
services_pols<-getting_models('perc_service', subset(pdf, type=='politician'))


models<-list(foreign_policy_pols, services_pols, protest_pols,
             criticism_pols)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'pols', 'Including Only Politicians')


#############################################
############Figure B7: Politicians Elected after 2010 Only


foreign_policy_pols2010<-getting_models('perc_foreign_policy', subset(pdf, type=='politician'&year_elected>=2010))
protest_pols2010<-getting_models('perc_protest', subset(pdf, type=='politician'&year_elected>=2010))
criticism_pols2010<-getting_models('perc_harsh_criticism', subset(pdf, type=='politician'&year_elected>=2010))
services_pols2010<-getting_models('perc_service', subset(pdf, type=='politician'&year_elected>=2010))

models<-list(foreign_policy_pols2010, services_pols2010,
             protest_pols2010, criticism_pols2010)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'pols2010', 'Including Only Politicians Elected after 2010')


############################
#####Figures B8 B9: Excluding/Controlling for Political Imprisonment


foreign_policy_expris<-getting_models('perc_foreign_policy', subset(pdf, prison==0))
protest_expris<-getting_models('perc_protest', subset(pdf, prison==0))
criticism_expris<-getting_models('perc_harsh_criticism', subset(pdf, prison==0))
services_expris<-getting_models('perc_service', subset(pdf, prison==0))

models<-list(foreign_policy_expris, services_expris, protest_expris, 
             criticism_expris)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'expris', 'Excluding Political Prisoners')



foreign_policy_conpris<-getting_models_polpris('perc_foreign_policy', pdf)
protest_conpris<-getting_models_polpris('perc_protest', pdf)
criticism_conpris<-getting_models_polpris('perc_harsh_criticism', pdf)
services_conpris<-getting_models_polpris('perc_service', pdf)

models<-list(foreign_policy_conpris, services_conpris,
             protest_conpris, criticism_conpris)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'conpris', 'Controlling for Political Imprisonment')



#################################
#######Figure B10: Controlling for party

foreign_policy_wparty<-getting_models_wparty('perc_foreign_policy', pdf)
protest_wparty<-getting_models_wparty('perc_protest', pdf)
criticism_wparty<-getting_models_wparty('perc_harsh_criticism', pdf)
services_wparty<-getting_models_wparty('perc_service', pdf)

models<-list(foreign_policy_wparty, services_wparty,
             protest_wparty, criticism_wparty)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'wparty', 'Controlling for Party')


#################################
#######Figure B11: Only VP

foreign_policy_onlyvp<-getting_models('perc_foreign_policy', subset(pdf, party=='VP'))
protest_onlyvp<-getting_models('perc_protest', subset(pdf, party=='VP'))
criticism_onlyvp<-getting_models('perc_harsh_criticism', subset(pdf, party=='VP'))
services_onlyvp<-getting_models('perc_service', subset(pdf, party=='VP'))


list_of_m<-list(foreign_policy_onlyvp, services_onlyvp,
                protest_onlyvp, criticism_onlyvp)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(list_of_m, dv_names, dv_types, model_types, 'onlyvp', 'Only Voluntad Popular Members')



################################
######Figs B12-B15: Leave one out results


####Leave-one-out test for coefficients

##Foreign policy

fp_df<-data.frame(matrix(ncol = 7, nrow = 0))
colnames(fp_df)<-c('term', 'estimate', 'std.error', 'statistic',
                   'p.value', 'conf.low', 'conf.high')

for (person in unique(pdf$actor.id)){
  es_foreign_policy <- plm(perc_foreign_policy ~ tweeted_exile+month+log(num_tweets+1), data = pdf, model = "within", subset=actor.id!=person)
  es_foreign_policy_coef <- coeftest(es_foreign_policy, vcov = vcovHC(es_foreign_policy, type = "HC1", cluster = "group"))
  coef_df<-subset(data.frame(tidy(es_foreign_policy_coef, coef.int=T)),term=='tweeted_exile') 
  coef_df$conf.low<-coef_df$estimate-coef_df$std.error*1.96
  coef_df$conf.high<-coef_df$estimate+coef_df$std.error*1.96
  
  fp_df<-rbind(fp_df, coef_df)
  
}

fp_df$id<-1:length(fp_df$estimate)

fp_loo<-ggplot(data=fp_df, aes(x=id, y=estimate))+
  geom_line()+ylim(0, 4.5)+
  geom_ribbon(data=fp_df, 
              aes(ymin=conf.low,ymax=conf.high), fill="gray", alpha=0.5)+
  theme_bw()+  theme(legend.position='bottom', axis.title=element_text(size=16),
                     text=element_text(size=16))+  xlab('Dropped User')+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  ylab('Change in Foreign Policy')+
  theme(plot.title = element_text(hjust = 0.5, size=18))

ggsave("plots/foreign_policy_loo.pdf", fp_loo, width = 8.5, height = 6.5)



##protest
pr_df<-data.frame(matrix(ncol = 7, nrow = 0))
colnames(pr_df)<-c('term', 'estimate', 'std.error', 'statistic',
                   'p.value', 'conf.low', 'conf.high')

for (person in unique(pdf$actor.id)){
  es_protest <- plm(perc_protest ~ tweeted_exile+month+log(num_tweets+1), data = pdf, model = "within", subset=actor.id!=person)
  es_protest_coef <- coeftest(es_protest, vcov = vcovHC(es_protest, type = "HC1", cluster = "group"))
  coef_df<-subset(data.frame(tidy(es_protest_coef, coef.int=T)),term=='tweeted_exile') 
  coef_df$conf.low<-coef_df$estimate-coef_df$std.error*1.96
  coef_df$conf.high<-coef_df$estimate+coef_df$std.error*1.96
  
  pr_df<-rbind(pr_df, coef_df)
  
}

pr_df$id<-1:length(pr_df$estimate)


pr_loo<-ggplot(data=pr_df, aes(x=id, y=estimate))+
  geom_line()+ylim(-1.5, 0)+
  geom_ribbon(data=pr_df, 
              aes(ymin=conf.low,ymax=conf.high), fill="gray", alpha=0.5)+
  theme_bw()+  theme(legend.position='bottom', axis.title=element_text(size=16),
                     text=element_text(size=16))+  xlab('Dropped User')+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  ylab('Change in Protest Discussion')+
  theme(plot.title = element_text(hjust = 0.5, size=18))

ggsave("plots/protest_loo.pdf", pr_loo, width = 8.5, height = 6.5)


#services
sv_df<-data.frame(matrix(ncol = 7, nrow = 0))
colnames(sv_df)<-c('term', 'estimate', 'std.error', 'statistic',
                   'p.value', 'conf.low', 'conf.high')

for (person in unique(pdf$actor.id)){
  es_sv <- plm(perc_service ~ tweeted_exile+month+log(num_tweets+1), data = pdf, model = "within", subset=actor.id!=person)
  es_sv_coef <- coeftest(es_sv, vcov = vcovHC(es_sv, type = "HC1", cluster = "group"))
  coef_df<-subset(data.frame(tidy(es_sv_coef, coef.int=T)),term=='tweeted_exile') 
  coef_df$conf.low<-coef_df$estimate-coef_df$std.error*1.96
  coef_df$conf.high<-coef_df$estimate+coef_df$std.error*1.96
  
  sv_df<-rbind(sv_df, coef_df)
  
}

sv_df$id<-1:length(sv_df$estimate)

sv_loo<-ggplot(data=sv_df, aes(x=id, y=estimate))+
  geom_line()+ylim(-4, 0)+
  geom_ribbon(data=sv_df, 
              aes(ymin=conf.low,ymax=conf.high), fill="gray", alpha=0.5)+
  theme_bw()+  theme(legend.position='bottom', axis.title=element_text(size=16),
                     text=element_text(size=16))+  xlab('Dropped User')+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  ylab('Change in Services Discussion')+
  theme(plot.title = element_text(hjust = 0.5, size=18))

ggsave("plots/services_loo.pdf", sv_loo, width = 8.5, height = 6.5)


#harsh criticism
hc_df<-data.frame(matrix(ncol = 7, nrow = 0))
colnames(hc_df)<-c('term', 'estimate', 'std.error', 'statistic',
                   'p.value', 'conf.low', 'conf.high')

for (person in unique(pdf$actor.id)){
  es_hc <- plm(perc_harsh_criticism ~ tweeted_exile+month+log(num_tweets+1), data = pdf, model = "within", subset=actor.id!=person)
  es_hc_coef <- coeftest(es_hc, vcov = vcovHC(es_hc, type = "HC1", cluster = "group"))
  coef_df<-subset(data.frame(tidy(es_hc_coef, coef.int=T)),term=='tweeted_exile') 
  coef_df$conf.low<-coef_df$estimate-coef_df$std.error*1.96
  coef_df$conf.high<-coef_df$estimate+coef_df$std.error*1.96
  
  hc_df<-rbind(hc_df, coef_df)
  
}

hc_df$id<-1:length(hc_df$estimate)

hc_loo<-ggplot(data=hc_df, aes(x=id, y=estimate))+
  geom_line()+ylim(0, 10)+
  geom_ribbon(data=hc_df, 
              aes(ymin=conf.low,ymax=conf.high), fill="gray", alpha=0.5)+
  theme_bw()+  theme(legend.position='bottom', axis.title=element_text(size=16),
                     text=element_text(size=16))+  xlab('Dropped User')+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  ylab('Change in Harsh Criticism')+
  theme(plot.title = element_text(hjust = 0.5, size=18))

ggsave("plots/criticism_loo.pdf", hc_loo, width = 8.5, height = 6.5)


##################################
#######Figure B16/B17: Goodman-Bacon decomposition

#Set up Data for Analysis
bacon_data<-balanced_panel
bacon_data$post<-bacon_data$tweeted_exile
bacon_data$year<-year(bacon_data$month)

#FOREIGN POLICY

#Group by actor.id (take mean for each year)
bacon_data_agg<-aggregate(x = bacon_data$perc_foreign_policy,                
                          by = list(bacon_data$actor.id, bacon_data$year, bacon_data$post),              
                          FUN = mean)
bacon_data_agg$foreign_policy_avg<-bacon_data_agg$x
bacon_data_agg$actor.id<-bacon_data_agg$Group.1
bacon_data_agg$year<-bacon_data_agg$Group.2
bacon_data_agg$post<-bacon_data_agg$Group.3
bacon_data_agg<-bacon_data_agg[c("post", "actor.id", "foreign_policy_avg", "year")]

#Remove Duplicates (pre & post for same year issue)
bacon_data_agg<-bacon_data_agg %>%
  distinct(year, actor.id, .keep_all = TRUE)

#Exclude people for which 2013 isn't post, otherwise
#"Treatment not weakly increasing with time" 
bacon_data_agg$zero2013<-ifelse(bacon_data_agg$year=="2013" & bacon_data_agg$post==1, 1,0)
table(bacon_data_agg$zero2013)
users_to_drop<-subset(bacon_data_agg, zero2013==1)
bacon_data_agg$drop<-ifelse(bacon_data_agg$actor.id %in% users_to_drop$actor.id, 1,0)
bacon_data_agg<-subset(bacon_data_agg, drop==0)

df_bacon <- bacon(foreign_policy_avg ~ post,
                  data = bacon_data_agg,
                  id_var = "actor.id",
                  time_var = "year")

coef_bacon <- sum(df_bacon$estimate * df_bacon$weight)
print(paste("Weighted sum of decomposition =", round(coef_bacon, 4)))

fit_tw <- lm(foreign_policy_avg ~ post + actor.id + factor(year), 
             data = bacon_data_agg)
print(paste("Two-way FE estimate =", round(fit_tw$coefficients[2], 4)))

ggplot(df_bacon) +
  aes(x = weight, y = estimate, shape = factor(type)) +
  labs(x = "Weight", y = "Estimate", shape = "Type") +
  geom_point()
ggsave("plots/bacon_goodman_foreign_policy.pdf", width = 11, height = 8)

#Make Table
webshot::install_phantomjs(force=TRUE)

df_bacon %>% 
  group_by(type) %>% 
  dplyr::summarize(Avg_Estimate = mean(estimate),
                   Number_Comparisons = n(),
                   Total_Weight = sum(weight)) %>% 
  kable(format = "html", booktabs = T, digits = 2, align = 'c',
        col.names = c("Type", "Average Estimate", "Number of 2x2 Comparisons", "Total Weight")) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"))%>%
  save_kable("plots/bacon_foreign_policy_table.pdf")

#PROTEST

#Group by actor.id (take mean for each year)
bacon_data_agg<-aggregate(x = bacon_data$perc_protest,                
                          by = list(bacon_data$actor.id, bacon_data$year, bacon_data$post),              
                          FUN = mean)
bacon_data_agg$protest_avg<-bacon_data_agg$x
bacon_data_agg$actor.id<-bacon_data_agg$Group.1
bacon_data_agg$year<-bacon_data_agg$Group.2
bacon_data_agg$post<-bacon_data_agg$Group.3
bacon_data_agg<-bacon_data_agg[c("post", "actor.id", "protest_avg", "year")]

#Remove Duplicates (pre & post for same year issue)
bacon_data_agg<-bacon_data_agg %>%
  distinct(year, actor.id, .keep_all = TRUE)

#Exclude people for which 2013 isn't post, otherwise
#"Treatment not weakly increasing with time" 
bacon_data_agg$zero2013<-ifelse(bacon_data_agg$year=="2013" & bacon_data_agg$post==1, 1,0)
table(bacon_data_agg$zero2013)
users_to_drop<-subset(bacon_data_agg, zero2013==1)
bacon_data_agg$drop<-ifelse(bacon_data_agg$actor.id %in% users_to_drop$actor.id, 1,0)
bacon_data_agg<-subset(bacon_data_agg, drop==0)

df_bacon <- bacon(protest_avg ~ post,
                  data = bacon_data_agg,
                  id_var = "actor.id",
                  time_var = "year")

coef_bacon <- sum(df_bacon$estimate * df_bacon$weight)
print(paste("Weighted sum of decomposition =", round(coef_bacon, 4)))

fit_tw <- lm(protest_avg ~ post + actor.id + factor(year), 
             data = bacon_data_agg)
print(paste("Two-way FE estimate =", round(fit_tw$coefficients[2], 4)))

ggplot(df_bacon) +
  aes(x = weight, y = estimate, shape = factor(type)) +
  labs(x = "Weight", y = "Estimate", shape = "Type") +
  geom_point()
ggsave("plots/bacon_goodman_protest.pdf", width = 11, height = 8)

#Make Table
df_bacon %>% 
  group_by(type) %>% 
  dplyr::summarize(Avg_Estimate = mean(estimate),
                   Number_Comparisons = n(),
                   Total_Weight = sum(weight)) %>% 
  kable(format = "html", booktabs = T, digits = 2, align = 'c',
        col.names = c("Type", "Average Estimate", "Number of 2x2 Comparisons", "Total Weight")) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"))%>%
  save_kable("plots/bacon_protest_table.pdf") 

#SERVICE PROVISION

#Group by actor.id (take mean for each year)
bacon_data_agg<-aggregate(x = bacon_data$perc_service,                
                          by = list(bacon_data$actor.id, bacon_data$year, bacon_data$post),              
                          FUN = mean)
bacon_data_agg$service_avg<-bacon_data_agg$x
bacon_data_agg$actor.id<-bacon_data_agg$Group.1
bacon_data_agg$year<-bacon_data_agg$Group.2
bacon_data_agg$post<-bacon_data_agg$Group.3
bacon_data_agg<-bacon_data_agg[c("post", "actor.id", "service_avg", "year")]

#Remove Duplicates (pre & post for same year issue)
bacon_data_agg<-bacon_data_agg %>%
  distinct(year, actor.id, .keep_all = TRUE)

#Exclude people for which 2013 isn't post, otherwise
#"Treatment not weakly increasing with time" 
bacon_data_agg$zero2013<-ifelse(bacon_data_agg$year=="2013" & bacon_data_agg$post==1, 1,0)
table(bacon_data_agg$zero2013)
users_to_drop<-subset(bacon_data_agg, zero2013==1)
bacon_data_agg$drop<-ifelse(bacon_data_agg$actor.id %in% users_to_drop$actor.id, 1,0)
bacon_data_agg<-subset(bacon_data_agg, drop==0)
bacon_data_agg$service_avg[is.na(bacon_data_agg$service_avg)]<-0

df_bacon <- bacon(service_avg ~ post,
                  data = bacon_data_agg,
                  id_var = "actor.id",
                  time_var = "year")

coef_bacon <- sum(df_bacon$estimate * df_bacon$weight)
print(paste("Weighted sum of decomposition =", round(coef_bacon, 4)))

fit_tw <- lm(service_avg ~ post + actor.id + factor(year), 
             data = bacon_data_agg)
print(paste("Two-way FE estimate =", round(fit_tw$coefficients[2], 4)))

ggplot(df_bacon) +
  aes(x = weight, y = estimate, shape = factor(type)) +
  labs(x = "Weight", y = "Estimate", shape = "Type") +
  geom_point()
ggsave("plots/bacon_goodman_service.pdf", width = 11, height = 8)

#Make Table
df_bacon %>% 
  group_by(type) %>% 
  dplyr::summarize(Avg_Estimate = mean(estimate),
                   Number_Comparisons = n(),
                   Total_Weight = sum(weight)) %>% 
  kable(format = "html", booktabs = T, digits = 4, align = 'c',
        col.names = c("Type", "Average Estimate", "Number of 2x2 Comparisons", "Total Weight")) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"))%>%
  save_kable("plots/bacon_service_table.pdf") 

#CRITICISM

#Group by actor.id (take mean for each year)
bacon_data_agg<-aggregate(x = bacon_data$perc_harsh_criticism,                
                          by = list(bacon_data$actor.id, bacon_data$year, bacon_data$post),              
                          FUN = mean)
bacon_data_agg$criticism_avg<-bacon_data_agg$x
bacon_data_agg$actor.id<-bacon_data_agg$Group.1
bacon_data_agg$year<-bacon_data_agg$Group.2
bacon_data_agg$post<-bacon_data_agg$Group.3
bacon_data_agg<-bacon_data_agg[c("post", "actor.id", "criticism_avg", "year")]

#Remove Duplicates (pre & post for same year issue)
bacon_data_agg<-bacon_data_agg %>%
  distinct(year, actor.id, .keep_all = TRUE)

#Exclude people for which 2013 isn't post, otherwise
#"Treatment not weakly increasing with time" 
bacon_data_agg$zero2013<-ifelse(bacon_data_agg$year=="2013" & bacon_data_agg$post==1, 1,0)
table(bacon_data_agg$zero2013)
users_to_drop<-subset(bacon_data_agg, zero2013==1)
bacon_data_agg$drop<-ifelse(bacon_data_agg$actor.id %in% users_to_drop$actor.id, 1,0)
bacon_data_agg<-subset(bacon_data_agg, drop==0)

df_bacon <- bacon(criticism_avg ~ post,
                  data = bacon_data_agg,
                  id_var = "actor.id",
                  time_var = "year")

coef_bacon <- sum(df_bacon$estimate * df_bacon$weight)
print(paste("Weighted sum of decomposition =", round(coef_bacon, 4)))

fit_tw <- lm(criticism_avg ~ post + actor.id + factor(year), 
             data = bacon_data_agg)
print(paste("Two-way FE estimate =", round(fit_tw$coefficients[2], 4)))

ggplot(df_bacon) +
  aes(x = weight, y = estimate, shape = factor(type)) +
  labs(x = "Weight", y = "Estimate", shape = "Type") +
  geom_point()
ggsave("plots/bacon_goodman_criticism.pdf", width = 11, height = 8)

#Make Table
df_bacon %>% 
  group_by(type) %>% 
  dplyr::summarize(Avg_Estimate = mean(estimate),
                   Number_Comparisons = n(),
                   Total_Weight = sum(weight)) %>% 
  kable(format = "html", booktabs = T, digits = 2, align = 'c',
        col.names = c("Type", "Average Estimate", "Number of 2x2 Comparisons", "Total Weight")) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"))%>%
  save_kable("plots/bacon_criticism_table.pdf") 


################################
#######Figures B18-B23: Callaway Sant'Anna Estimator

data_panel$first.treat<-(interval(rep(as.Date('2013-01-01'), length(data_panel$month_of_exile)) ,data_panel$month_of_exile) %/% months(1))
data_panel$first.treat[is.na(data_panel$month_of_exile)]<-90
data_panel$month_num<-(interval(rep(as.Date('2013-01-01'), length(data_panel$month)) ,data_panel$month) %/% months(1))
data_panel$id<-as.numeric(as.factor(data_panel$actor.id))

fp_att <- att_gt(yname = "perc_foreign_policy",
                 tname = "month_num",
                 idname = "id",
                 gname = "first.treat",
                 allow_unbalanced_panel = TRUE,
                 xformla = ~1,
                 data = data_panel)

fp_att_dyn <- aggte(fp_att, type = "dynamic", na.rm = TRUE)
summary(fp_att_dyn)

ggsave('plots/ggdid_fp.pdf', ggdid(fp_att_dyn, xgap=10)+theme(axis.title=element_text(size=22),
                                                              axis.text=element_text(size=22),
                                                              plot.title=element_text(size=22),
                                                              legend.text=element_text(size=22)), width=9, height=7)

att<-fp_att_dyn$overall.att
seatt<-fp_att_dyn$overall.se
clow<-att-1.96*seatt
chigh<-att+1.96*seatt

calsant<-data.frame(cbind('Foreign Policy', att, seatt, clow, chigh))

sv_att <- att_gt(yname = "perc_service",
                 tname = "month_num",
                 idname = "id",
                 gname = "first.treat",
                 allow_unbalanced_panel = TRUE,
                 xformla = ~1,
                 data = data_panel)


sv_att_dyn <- aggte(sv_att, type = "dynamic", na.rm = TRUE)

summary(sv_att_dyn)

att<-sv_att_dyn$overall.att
seatt<-sv_att_dyn$overall.se
clow<-att-1.96*seatt
chigh<-att+1.96*seatt


ggsave('plots/ggdid_sv.pdf', ggdid(sv_att_dyn, xgap=10)+theme(axis.title=element_text(size=22),
                                                              axis.text=element_text(size=22),
                                                              plot.title=element_text(size=22),
                                                              legend.text=element_text(size=22)), width=9, height=7)

calsant<-rbind(calsant, data.frame(cbind('Service Provision', att, seatt, clow, chigh)))


pr_att <- att_gt(yname = "perc_protest",
                 tname = "month_num",
                 idname = "id",
                 gname = "first.treat",
                 allow_unbalanced_panel = TRUE,
                 xformla = ~1,
                 data = data_panel)


pr_att_dyn <- aggte(pr_att, type = "dynamic", na.rm = TRUE)


att<-pr_att_dyn$overall.att
seatt<-pr_att_dyn$overall.se
clow<-att-1.96*seatt
chigh<-att+1.96*seatt

ggsave('plots/ggdid_pr.pdf', ggdid(pr_att_dyn, xgap=10)+theme(axis.title=element_text(size=22),
                                                              axis.text=element_text(size=22),
                                                              plot.title=element_text(size=22),
                                                              legend.text=element_text(size=22)), width=9, height=7)

calsant<-rbind(calsant, data.frame(cbind('Protest', att, seatt, clow, chigh)))


hc_att_all <- att_gt(yname = "perc_harsh_criticism",
                     tname = "month_num",
                     idname = "id",
                     gname = "first.treat",
                     allow_unbalanced_panel = TRUE,
                     xformla = ~1,
                     data = data_panel)


hc_att_dyn_all <- aggte(hc_att_all, type = "dynamic", na.rm = TRUE)

summary(hc_att_dyn_all)

att<-hc_att_dyn_all$overall.att
seatt<-hc_att_dyn_all$overall.se
clow<-att-1.96*seatt
chigh<-att+1.96*seatt


ggsave('plots/ggdid_hc_all.pdf', ggdid(hc_att_dyn_all, xgap=10)+theme(axis.title=element_text(size=22),
                                                                      axis.text=element_text(size=22),
                                                                      plot.title=element_text(size=22),
                                                                      legend.text=element_text(size=22)), width=9, height=7)

calsant<-rbind(calsant, data.frame(cbind('Harsh Criticism', att, seatt, clow, chigh)))


hc_att <- att_gt(yname = "perc_harsh_criticism",
                 tname = "month_num",
                 idname = "id",
                 gname = "first.treat",
                 allow_unbalanced_panel = TRUE,
                 xformla = ~1,
                 data = subset(data_panel, year(month)>2015))


hc_att_dyn <- aggte(hc_att, type = "dynamic", na.rm = TRUE)

summary(hc_att_dyn)
ggsave('plots/ggdid_hc_p2018.pdf', ggdid(hc_att_dyn, xgap=10)+theme(axis.title=element_text(size=22),
                                                                    axis.text=element_text(size=22),
                                                                    plot.title=element_text(size=22),
                                                                    legend.text=element_text(size=22)), width=9, height=7)

att<-hc_att_dyn$overall.att
seatt<-hc_att_dyn$overall.se
clow<-att-1.96*seatt
chigh<-att+1.96*seatt


calsant<-rbind(calsant, data.frame(cbind('Harsh Criticism 2016-', att, seatt, clow, chigh)))

level_order<-c('Harsh Criticism 2016-', 'Harsh Criticism', 'Protest', 'Service Provision', 'Foreign Policy')

callsant<-ggplot(calsant, aes(x=factor(V1, level = level_order), y=as.numeric(att)))+
  geom_pointrange(aes(ymin=as.numeric(clow), ymax=as.numeric(chigh)), position=position_dodge(width=.3))+
  theme_bw()+
  theme(legend.position='bottom', axis.title=element_text(size=16),
        text=element_text(size=18))+
  coord_flip()+
  geom_hline(aes(yintercept=0), lty=4)+
  scale_y_continuous(name='% change in tweets')+
  scale_shape_discrete(name=' ', guide=guide_legend(reverse=T))+
  scale_linetype_discrete(name=' ', guide=guide_legend(reverse=T))+
  xlab("")+
  theme(plot.title = element_text(hjust = 0.5, size=20),
        axis.title=element_text(size=22),
        axis.text=element_text(size=20))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
ggsave("plots/callsant.pdf", callsant, width = 7.5, height = 6.5)


######################################################
##############Figure B24: Excluding those exiled before 2013
foreign_policy_ep2013<-getting_models('perc_foreign_policy', subset(pdf, exile=='no'|date_of_exile>='2013-01-01'))
protest_ep2013<-getting_models('perc_protest', subset(pdf, exile=='no'|date_of_exile>='2013-01-01'))
criticism_ep2013<-getting_models('perc_harsh_criticism', subset(pdf, exile=='no'|date_of_exile>='2013-01-01'))
services_ep2013<-getting_models('perc_service', subset(pdf, exile=='no'|date_of_exile>='2013-01-01'))

models<-list(foreign_policy_ep2013,services_ep2013,
             protest_ep2013, criticism_ep2013)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'ep2013', 'Exiled After 2013')

#############################################################
#############Figure B25: Excluding those exiled before 2017
foreign_policy_ep2017<-getting_models('perc_foreign_policy', subset(pdf, exile=='no'|date_of_exile>='2017-01-01'))
protest_ep2017<-getting_models('perc_protest', subset(pdf, exile=='no'|date_of_exile>='2017-01-01'))
criticism_ep2017<-getting_models('perc_harsh_criticism', subset(pdf, exile=='no'|date_of_exile>='2017-01-01'))
services_ep2017<-getting_models('perc_service', subset(pdf, exile=='no'|date_of_exile>='2017-01-01'))

models<-list(foreign_policy_ep2017, services_ep2017, 
             protest_ep2017, criticism_ep2017)

dv_names<-c('Foreign Policy', 'Services', 'Protest', 'Criticism')

dv_types<-c('bold', 'bold', 'bold', 'bold')

model_types<-c('Two-Way FEs', 'Month FEs', 'Time Trend')

coef_plot(models, dv_names, dv_types, model_types, 'ep2017', 'Exiled After 2017')


